## ---------------------------
##
## Script name: MRW_high_hurdles_main.R
##
## Purpose of script: Replicate main findings from "High Hurdles: Legislative 
##                    Professionalism and the Effectiveness of Women State 
##                    Legislators" - 2025, JOP
##
## Author: Rob McGrath, Josh Ryan, Tia Wrighten
##
## Date Created: 2025-04-08
##
## ---------------------------


rm(list=ls())

if (!require("pacman")) install.packages("pacman")
pacman::p_load(rstudioapi, tidyverse, plm, AER, texreg, broom, readr)

my.dir <- dirname(rstudioapi::getActiveDocumentContext()$path)
setwd(my.dir)

d <- read_csv("legislators.csv")


######### TABLE 1 ##############

d %>% 
  plm(log_SLES ~ female + year, data=., 
      index=c("districtFE"), model="within") -> sles_basic_dfe1

d %>% 
  plm(log_SLES ~ female  + factor(year), data=., 
      index=c("districtFE"), model="within") -> sles_basic_dfe2

d %>% 
  plm(log_SLES ~ female   + factor(st)*year, data=., 
      index=c("districtFE"), model="within") -> sles_basic_dfe3

d %>%   filter(st!= "CT" & st!="MA" & st!="KY" & st!="NE") %>% 
  plm(log_SLES ~ female*prop_fem_chamber + female*female_leader + 
        female*log_rawcommvalues_scaled + 
        female*ncsl  +
        yearsserved + mp_member  + ideological_deviance  + year, data=., 
      index=c("districtFE"), model="within") -> sles_basic_dfe4

d %>%    filter(st!= "CT" & st!="MA" & st!="KY" & st!="NE") %>% 
  plm(log_SLES ~ female*prop_fem_chamber + female*female_leader + 
        female*log_rawcommvalues_scaled +
        female*ncsl +
        yearsserved + mp_member  + ideological_deviance  +
        factor(year), data=., 
      index=c("districtFE"), model="within") -> sles_basic_dfe5

d %>%   filter(st!= "CT" & st!="MA" & st!="KY" & st!="NE") %>% 
  plm(log_SLES ~ female*prop_fem_chamber + female*female_leader + 
        female*log_rawcommvalues_scaled +
        female*ncsl  +
        yearsserved + mp_member  + ideological_deviance  + 
        factor(st)*year, data=., 
      index=c("districtFE"), model="within") -> sles_basic_dfe6

fits <- list(sles_basic_dfe1, sles_basic_dfe2, sles_basic_dfe3, sles_basic_dfe4, sles_basic_dfe5, sles_basic_dfe6)

coefs = list(
  "female" = "Female",
  "prop_fem_chamber" = "Female Proportion in Chamber (Scaled to Mean 0)",
  "female_leader" = "Female Leader",
  "log_rawcommvalues_scaled" = "Committee Portfolio Value (Logged) - Scaled",
  "female:prop_fem_chamber" = "Female x Scaled Female Proportion in Chamber",
  "female:female_leader" = "Female x Female Leader",
  "female:ncslGray" = "Female x Hybrid Legislature",
  "female:ncslGreen" = "Female x Professional Legislature",
  "female:log_rawcommvalues_scaled" = "Female x Scaled Committee Portfolio Value",
  "yearsserved" = "Seniority (in Years)",
  "mp_member" = "Majority Party Member",
  "democrat" = "Member is a Democrat",
  "Leader" = "Leader",
  "ideological_deviance" = "Ideological Distance from Party Median",
  "(Intercept)" = "Intercept" 
)

texreg(list(coeftest(sles_basic_dfe1, vcov=vcovHC(sles_basic_dfe1,type="HC0",cluster="group")),
            coeftest(sles_basic_dfe2, vcov=vcovHC(sles_basic_dfe2,type="HC0",cluster="group")),
            coeftest(sles_basic_dfe3, vcov=vcovHC(sles_basic_dfe3,type="HC0",cluster="group")),
            coeftest(sles_basic_dfe4, vcov=vcovHC(sles_basic_dfe4,type="HC0",cluster="group")),
            coeftest(sles_basic_dfe5, vcov=vcovHC(sles_basic_dfe5,type="HC0",cluster="group")),
            coeftest(sles_basic_dfe6, vcov=vcovHC(sles_basic_dfe6,type="HC0",cluster="group"))),
       custom.coef.map = coefs,
       custom.gof.rows = list(
         "N" = c(sapply(fits,function(x) glance(x)$nobs)),
         "District Fixed Effects" = c("yes", "yes", "yes", "yes", "yes", "yes"),
         "Year Fixed Effects" = c("no", "yes", "no", "no", "yes", "no"),
         "State x Year trend" = c("no", "no", "yes", "no", "no", "yes")#, 
       ),
       custom.model.names = c("1", "2", "3", "4", "5", "6"),
       caption = "Sex-Based Differences in Logged State Legislative Effectiveness Scores, 2007-2018",
       table=TRUE,
       digits = 3)




########### FIGURE 1 ################


model <- sles_basic_dfe5

cat("Marginal Effect of Female on Legislative Effectiveness for Amateur Legislatures: ", model$coef["female"])
cat("Marginal Effect of Female on Legislative Effectiveness for Hybrid Legislatures: ", model$coef["female"] + model$coef["female:ncslGray"] )
cat("Marginal Effect of Female on Legislative Effectiveness for Professional Legislatures: ", model$coef["female"] +
      model$coef["female:ncslGreen"] )

coef.vec <- c(model$coef["female"], model$coef["female"] + model$coef["female:ncslGray"], 
              model$coef["female"] + model$coef["female:ncslGreen"])

model.vcov <- vcov(model)

se.vec <- c(sqrt(model.vcov["female","female"]), 
            sqrt(model.vcov["female","female"] + model.vcov["female:ncslGray","female:ncslGray"] +
                   2*model.vcov["female","female:ncslGray"]),
            sqrt(model.vcov["female","female"] + model.vcov["female:ncslGreen","female:ncslGreen"] + 
                   2*model.vcov["female","female:ncslGreen"]))


# Extract robust variance-covariance matrix
robust_vcov <- vcovHC(model, type = "HC0", cluster = "group")

# Compute robust standard errors for the coefficients of interest
robust_se <- c(
  sqrt(robust_vcov["female", "female"]),
  sqrt(
    robust_vcov["female", "female"] +
      robust_vcov["female:ncslGray", "female:ncslGray"] +
      2 * robust_vcov["female", "female:ncslGray"]
  ),
  sqrt(
    robust_vcov["female", "female"] +
      robust_vcov["female:ncslGreen", "female:ncslGreen"] +
      2 * robust_vcov["female", "female:ncslGreen"]
  )
)



# Create a data frame for plotting with reordered levels
plot_data <- data.frame(
  Legislature = factor(
    c("Amateur Legislatures", "Hybrid Legislatures", "Professional Legislatures"),
    levels = c("Professional Legislatures", "Hybrid Legislatures", "Amateur Legislatures") # Change the order here
  ),
  Coefficient = coef.vec,
  Lower_CI = coef.vec - qnorm(0.95) * robust_se,
  Upper_CI = coef.vec + qnorm(0.95) * robust_se
)

# Create the plot
ggplot(plot_data, aes(x = Coefficient, y = Legislature)) +
  geom_point(size = 3) +  # Add points for the coefficients
  geom_errorbarh(aes(xmin = Lower_CI, xmax = Upper_CI), height = 0.2) +  # Add horizontal error bars
  geom_vline(xintercept = 0, linetype = "dashed") +  # Add a vertical reference line at 0
  labs(
    title = "Marginal Effects of Female, with 95% CI",
    x = "Marginal Effect",
    y = ""
  ) +
  theme_minimal() +  # Use a clean theme
  theme(
    axis.text.y = element_text(size = 10),
    axis.title.x = element_text(size = 12),
    plot.title = element_text(hjust = 0.5, size = 14)
  )




######### TABLE 2 ##############

d %>%   filter(
  st != "CT" & st != "MA" & st != "KY" & st != "NE", 
  !is.na(log_SLES)                                  
) %>%
  plm(log_rawcommvalues ~ female + year, data=., 
      index=c("districtFE"), model="within") -> pv_basic_dfe1

d %>%   filter(
  st != "CT" & st != "MA" & st != "KY" & st != "NE", 
  !is.na(log_SLES)                                  
) %>%
  plm(log_rawcommvalues ~ female + factor(year), data=., 
      index=c("districtFE"), model="within") -> pv_basic_dfe2

d %>% filter(
  st != "CT" & st != "MA" & st != "KY" & st != "NE", 
  !is.na(log_SLES)                                  
) %>%
  plm(log_rawcommvalues ~ female + factor(st)*year, data=., 
      index=c("districtFE"), model="within") -> pv_basic_dfe3

d %>% filter(
  st != "CT" & st != "MA" & st != "KY" & st != "NE", 
  !is.na(log_SLES)                                  
) %>%
  plm(log_rawcommvalues ~ female*prop_fem_chamber + female*female_leader + female*ncsl + 
        yearsserved + mp_member + ideological_deviance +
        year, data=., 
      index=c("districtFE"), model="within") -> pv_basic_dfe4

d %>% filter(
  st != "CT" & st != "MA" & st != "KY" & st != "NE", 
  !is.na(log_SLES)                                  
) %>%
  plm(log_rawcommvalues ~ female*prop_fem_chamber + female*female_leader + female*ncsl + 
        yearsserved + mp_member + ideological_deviance +
        factor(year), data=., 
      index=c("districtFE"), model="within") -> pv_basic_dfe5

d %>% filter(
  st != "CT" & st != "MA" & st != "KY" & st != "NE", 
  !is.na(log_SLES)                                  
) %>%
  plm(log_rawcommvalues ~ female*prop_fem_chamber + female*female_leader + female*ncsl + 
        yearsserved + mp_member + ideological_deviance +
        factor(st)*year, data=., 
      index=c("districtFE"), model="within") -> pv_basic_dfe6

fits <- list(pv_basic_dfe1, pv_basic_dfe2, pv_basic_dfe3, pv_basic_dfe4, pv_basic_dfe5, pv_basic_dfe6)

coefs = list(
  "female" = "Female",
  "prop_fem_chamber" = "Female Proportion in Chamber (Scaled to Mean 0)",
  "female_leader" = "Female Leader",
  "female:prop_fem_chamber" = "Female x Scaled Female Proportion in Chamber",
  "female:female_leader" = "Female x Female Leader",
  "female:ncslGray" = "Female x Hybrid Legislature",
  "female:ncslGreen" = "Female x Professional Legislature",
  "yearsserved" = "Seniority (in Years)",
  "mp_member" = "Majority Party Member",
  "ideological_deviance" = "Ideological Distance from Party Median",
  "(Intercept)" = "Intercept" 
)

texreg(list(coeftest(pv_basic_dfe1, vcov=vcovHC(pv_basic_dfe1,type="HC0",cluster="group")),
            coeftest(pv_basic_dfe2, vcov=vcovHC(pv_basic_dfe2,type="HC0",cluster="group")),
            coeftest(pv_basic_dfe3, vcov=vcovHC(pv_basic_dfe3,type="HC0",cluster="group")),
            coeftest(pv_basic_dfe4, vcov=vcovHC(pv_basic_dfe4,type="HC0",cluster="group")),
            coeftest(pv_basic_dfe5, vcov=vcovHC(pv_basic_dfe5,type="HC0",cluster="group")),
            coeftest(pv_basic_dfe6, vcov=vcovHC(pv_basic_dfe6,type="HC0",cluster="group"))),
       custom.coef.map = coefs,
       custom.gof.rows = list(
         "N" = c(sapply(fits,function(x) glance(x)$nobs)),
         "District Fixed Effects" = c("yes", "yes", "yes", "yes", "yes", "yes"),
         "Year Fixed Effects" = c("no", "yes", "no", "no", "yes", "no"),
         "State x Year trend" = c("no", "no", "yes", "no", "no", "yes")#, 
       ),
       custom.model.names = c("1", "2", "3", "4", "5", "6"),
       caption = "Sex-Based Differences in Logged Committee Portfolio Values, 2007-2018",
       table=TRUE,
       digits = 3)





